Develop TAZ/Census Equivalency Table

Single TAZ to Tract Equivalence in New Tracts

Lets first identify all tracts in sets (a) and (b)

These tracts should have an area of intersection with a TAZ which is at least as large (within error) as the TAZ.

We say “within error” because the topology errors here make full identity unlikely.

Ideally, the Tract’s area of intersection with the TAZ would be identical to the TAZ’s area (or greater). But because of topological errors, this rarely occurs. So we have to decide what cutoff we want to say is “within error.”

Below we set it to conditions where the Tract’s area of intersection with the TAZ is at least 70% of the area of the TAZ.

We can compare this to the ratio of single-taz intersections (types (a) and (b)) in the year 2000 data.

intersection_df_new_tracts$intersection_area_over_tract_area <- as.numeric(intersection_df_new_tracts$intersection_area)/as.numeric(intersection_df_new_tracts$tract_area) 
summary(intersection_df_new_tracts$intersection_area_over_tract_area)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0001366 0.3553154 0.9610558 0.7221450 0.9862458 1.0000000

print(table(intersection_df_new_tracts$intersection_area_over_tract_area>.70))

## 
## FALSE  TRUE 
##   146   363

If we base our decision on what on 2000 data then we should expect about 96% of the Tracts to map to only 1 TAZ, which again would be the relationship conditions defined as (a) and (b).

At 70% of intersection we get 363/387 or 93% of the intersections in the new Tracts being single-taz intersections.

This seems about right, and a bit conservative.

Interestingly, setting the intersection_area_over_tract_area ratio to 50% gets us to almost exactly 96%, which is what the year 2000 data had. But lets be conservative and go with 70%

Lets call those intersections acceptable, and put them in their own data frame.

intersection_type_a_b <- intersection_df_new_tracts[intersection_df_new_tracts$intersection_area_over_tract_area>.70,]

tracts_left_id <- intersection_df_new_tracts[!intersection_df_new_tracts$tract %in% intersection_type_a_b$tract,]$tract
length(tracts_left_id)

## [1] 35

We’ll leave that data frame (intersection_type_a_b) aside and add them to the main table later with the type (c) intersections.

Multiple TAZ to Tract Equivalence in New Tracts

That leaves us just a handful of tracts left. We’ll put the spatial intersections for these types (c) into their own data frame.

intersection_tracts_left <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tracts_left_id,]
dim(intersection_tracts_left)

## [1] 35 15

rm(intersection_tracts_left)

Lets have a look at them:

tract_ids <- tracts_left_id

taz_overlap_ids <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]$taz

m_tracts10 <- tracts10[tracts10$tract %in% tract_ids,]
m_taz <- taz_df[taz_df$taz %in% taz_overlap_ids,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]

mapview(m_tracts10, col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
  mapview(m_taz, col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
  mapview(m_intersection, color="red", col.regions="red", alpha=0.8, map.types=c('Stamen.Toner.Light'))

Below we’ll review a few of these intersections in more detail in order to determine which we should keep.

TAZ larger than Tract

Others are examples of where the TAZ is larger than the Tract. For example:

cross_tract <- "505009"
tract_ids <- cross_tract

taz_overlap_ids <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]$taz

m_tracts10 <- tracts10[tracts10$tract %in% tract_ids,]
m_taz <- taz_df[taz_df$taz %in% taz_overlap_ids,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]

mapview(m_tracts10, col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
  mapview(m_taz, col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
  mapview(m_intersection, color="red", col.regions="red", alpha=0.8, map.types=c('Stamen.Toner.Light'))

This is a type (b) and type (c) coverage. Some of the remaining tracts are like this. From here on out we’ll just call the intersections valid or not.

Tract Is Split by TAZ

In some cases the Tract seems to be split by a TAZ.

cross_tract <- "403502"
tract_ids <- cross_tract
#m_blocks <- blocks[blocks$TRACTCE10 %in% tract_ids,]

taz_overlap_ids <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]$taz

m_tracts10 <- tracts10[tracts10$tract %in% tract_ids,]
m_taz <- taz_df[taz_df$taz %in% taz_overlap_ids,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_ids,]

mapview(m_tracts10, col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
  mapview(m_taz, col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
  mapview(m_intersection, color="red", col.regions="red", alpha=0.8, map.types=c('Stamen.Toner.Light'))

We might want to develop rules for this kind of TAZ/Tract relationship in the future. For now, we will just say its equivalent to the 1 TAZ that majority covers it (972).

intersection_type_a_b <- rbind(intersection_type_a_b,
                                intersection_df_new_tracts[intersection_df_new_tracts$tract=="403502" &
                                  intersection_df_new_tracts$taz==972,])

Manually Identify Valid Multi-Intersection Tracts

Valid intersections include those in tracts: 511705, 614000, 061500, 061100, 131100, 410500, 313206, 061200,601902, 601901

990100 is not valid, because it is a park/open space preserve. 980401 is not valid because its the Farallon islands.

Lets add the valid intersections to the table.

valid_intersections <- c("505009","613900","061500","511705", "614000", 
                         "061500", "061100", "131100", "410500", 
                         "313206","061200","601902", "601901")
intersection_type_c <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% valid_intersections,]

intersection_type_c <- intersection_type_c

st_geometry(intersection_type_c) <- NULL

intersection_type_c <- dplyr::select(intersection_type_c,taz,tract)

intersection_type_c <- intersection_type_c[,c('tract','taz')]

intersection_type_c$header_string <- 'rtaz'

intersection_type_c$num <- ave(intersection_type_c[['taz']], 
              intersection_type_c[['tract']], 
              FUN = seq_along)

intersection_type_c_sparse <- intersection_type_c %>% 
  tidyr::unite("header_string", 
                   header_string, 
                   num) %>% 
    tidyr::spread(header_string, taz)
#rm(intersection_type_c)

Adding a Column to the Lookup Table

When we build the sparse table for multiple intersections, it becomes apparent that we will have to add a column to the lookup table. In 2000, there were 5 columns, but here we’ve found that (for 2010) there are 6 valid intersections for 1 tract. That tract is in downtown SF, and a map of it can be found in the Appendix.

Lets take a look at what we’re adding to the updated 2000 table before we do it. We may need to reshape the columns.

st_geometry(intersection_type_a_b) <- NULL

intersection_type_a_b <- intersection_type_a_b[,c('tract','taz')]

knitr::kable(head(intersection_type_c_sparse))

knitr::kable(head(intersection_type_a_b))

knitr::kable(head(Tract_zone_2000))

Reshape (a), (b), (c) Equivalence then Add to with Year 2000

So we need to reshape the type_a_b table and add the new column to the 2000 table.

And rename columns to match (underscores, etc)

Tract_zone_2010 <- Tract_zone_2000
intersection_type_a_b <- dplyr::rename(intersection_type_a_b, rtaz1 = taz)
intersection_type_a_b$rtaz2 <- as.integer(NA)
intersection_type_a_b$rtaz3 <- as.integer(NA)
intersection_type_a_b$rtaz4 <- as.integer(NA)
intersection_type_a_b$rtaz5 <- as.integer(NA)
intersection_type_a_b$rtaz6 <- as.integer(NA)
Tract_zone_2010$rtaz6 <- as.integer(NA)
names(intersection_type_a_b) <- names(Tract_zone_2010)
names(intersection_type_c_sparse) <- names(Tract_zone_2010)

Now we can bind them all together. Lets have a look at them and then do so.

knitr::kable(head(Tract_zone_2010))

knitr::kable(head(intersection_type_a_b))

knitr::kable(head(intersection_type_c_sparse))

Tract_zone_2010 <- rbind(Tract_zone_2010,intersection_type_a_b,intersection_type_c_sparse)
dim(Tract_zone_2010)

## [1] 1781    7

#rm(intersection_type_a_b)
#rm(intersection_type_c_sparse)

We have too many tracts now. Thats because we still have tracts that were changed from 2000 to 2010 in the table that we have not dropped.

Drop the tracts updated for 2010 (by id)

Tract_zone_2010 <- Tract_zone_2010[!Tract_zone_2010$tract %in% tracts_to_drop,]
dim(Tract_zone_2010)

## [1] 1581    7

dim(tracts10)

## [1] 1592   10

We still have a few missing.

Map those tracts that are missing from the lookup.

missing_tracts10 <- tracts10[!tracts10$tract %in% Tract_zone_2010$tract,]

Check Missing Tracts are in Water

These should just be water tracts. We will load tracts not cut by water and water here to double check that.

options(tigris_use_cache=TRUE)
library(tigris)
counties1<-c("01","13","41","55","75","81","85","95","97")
tracts_not_cut <- tigris::tracts("CA", counties1, class="sf", year=2010)
detach("package:tigris", unload=TRUE)

tracts_not_cut <- dplyr::select(tracts_not_cut,TRACTCE10)
tracts_not_cut <- dplyr::rename(tracts_not_cut,tract = TRACTCE10)
tracts_not_cut <- st_transform(tracts_not_cut, crs=26910)
tracts_not_cut$tract_total_area <- st_area(tracts_not_cut)

bay_water <- st_read("https://geo.nyu.edu/download/file/stanford-mb777jk0330-geojson.json")

## Reading layer `stanford-mb777jk0330-geojson' from data source `https://geo.nyu.edu/download/file/stanford-mb777jk0330-geojson.json' using driver `GeoJSON'
## Simple feature collection with 367 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -123.6336 ymin: 36.92222 xmax: -121.2531 ymax: 38.85059
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

bay_water <- bay_water[st_is_valid(bay_water),]
bay_water <- st_transform(bay_water, crs=26910)

df1 <- st_intersection(tracts_not_cut, bay_water)

df1$intersection_area <- st_area(df1)

mapview(df1[df1$tract %in% missing_tracts10$tract,])

rm(bay_water)
rm(tracts_not_cut)
rm(df1)

Fill NA Values for Water Tracts in Table

These tracts are all over water.

We can put them into the table for clarity

water_tracts10 <- data.frame(tract=c(unique(as.character(missing_tracts10$tract))))
water_tracts10$rtaz1 <- as.integer(NA)
water_tracts10$rtaz2 <- as.integer(NA)
water_tracts10$rtaz3 <- as.integer(NA)
water_tracts10$rtaz4 <- as.integer(NA)
water_tracts10$rtaz5 <- as.integer(NA)
water_tracts10$rtaz6 <- as.integer(NA)
knitr::kable(head(water_tracts10))

Tract_zone_2010 <- rbind(Tract_zone_2010,water_tracts10)
dim(Tract_zone_2010)[1]==length(unique(tracts10$tract))

## [1] TRUE

rm(water_tracts10)

It has the right number of rows

Do all TAZ’s have a tract?

tract_zone_2010_dense <- reshape2::melt(Tract_zone_2010, id.vars = "tract", na.rm=TRUE)
tazs <- unique(tract_zone_2010_dense$value)
table(taz_df$taz %in% tazs)

## 
## FALSE  TRUE 
##     8  1446

Add Valid Intersections not in the 2000 Table

There are 8 TAZ’s not in the table. Were they in the intersection types a, b or c we identified?

missing_tazs <- taz_df[!taz_df$taz %in% tazs,]$taz
table(missing_tazs %in% intersection_type_a_b$rtaz1)

## 
## FALSE 
##     8

table(missing_tazs %in% intersection_type_c$taz)

## 
## FALSE 
##     8

They were not.

Are they not new tracts?

table(intersection_df[intersection_df$taz %in% missing_tazs,]$tract %in% tracts00$tract)

## 
## FALSE  TRUE 
##    11    13

They were not because they were in the year 2000 tract table.

tract_zone_2000_dense[tract_zone_2000_dense$taz %in% missing_tazs,]

## [1] tract    variable value   
## <0 rows> (or 0-length row.names)

They were not in the 2000 equivalency table.

Where are they?

tract_overlap_ids <- intersection_df[intersection_df$taz %in% missing_tazs,]$tract

m_tracts10 <- tracts10[tracts10$tract %in% tract_overlap_ids,]
m_taz <- taz_df[taz_df$taz %in% missing_tazs,]
m_intersection <- intersection_df_new_tracts[intersection_df_new_tracts$tract %in% tract_overlap_ids,]

mapview(m_tracts10, col.regions="green", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
  mapview(m_taz, col.regions="blue", alpha=0.8, map.types=c('Stamen.Toner.Light')) +
  mapview(m_intersection, color="red", col.regions="red", alpha=0.8, map.types=c('Stamen.Toner.Light'))

These seem like intersections that should be included based on our definitions. All of these would fall into both types (b) and (c) intersections.

So we’ll add them to the table

additions_to_2000_lookup <- m_intersection[,c('tract','taz')]
st_geometry(additions_to_2000_lookup) <- NULL

Tract_zone_2010_temp <- Tract_zone_2010[!Tract_zone_2010$tract %in% additions_to_2000_lookup$tract,]

additions_to_2000_lookup$header_string <- 'rtaz'

additions_to_2000_lookup$num <- ave(additions_to_2000_lookup[['taz']], 
              additions_to_2000_lookup[['tract']], 
              FUN = seq_along)

additions_to_2000_lookup_sparse <- additions_to_2000_lookup %>% 
  tidyr::unite("header_string", 
                   header_string,
                   sep="",
                   num) %>% 
    tidyr::spread(header_string, taz)

additions_to_2000_lookup_sparse$rtaz6 <- as.integer(NA)

Tract_zone_2010_temp <- rbind(Tract_zone_2010_temp,additions_to_2000_lookup_sparse)

dim(Tract_zone_2010_temp)

## [1] 1584    7

Final Checks

Check again for missing TAZ’s:

tract_zone_2010_dense <- reshape2::melt(Tract_zone_2010_temp, id.vars = "tract", na.rm=TRUE)
tazs <- unique(tract_zone_2010_dense$value)
table(taz_df$taz %in% tazs)

## 
## TRUE 
## 1454

No missing TAZ’s.

missing_tracts <- tracts10[!tracts10$tract %in% Tract_zone_2010$tract,]
dim(missing_tracts)

## [1]  0 10

No missing tracts.

TAZ:

mapview(m_taz, color="blue", col.regions="transparent", alpha=1, map.types=c('Stamen.Toner.Light'))

Intersection:

mapview(m_intersection,color="green", col.regions="transparent", alpha=1, map.types=c('Stamen.Toner.Light'))

Blocks:

mapview(m_blocks, color="magenta", col.regions="transparent", alpha=1, map.types=c('Stamen.Toner.Light'))

Tracts, TAZ, and Blocks:

mapview(m_tracts10, color="red", col.regions="transparent", alpha=0.75, map.types=c('Stamen.Toner.Light')) +
  mapview(m_blocks, color="magenta", col.regions="transparent", alpha=0.75, map.types=c('Stamen.Toner.Light')) +
  mapview(m_taz, color="blue", col.regions="transparent", alpha=0.75, map.types=c('Stamen.Toner.Light'))

Use Blocks, not Geometries?

If TAZ boundaries can be built from Census Blocks we might need not have to draw geometries at all.

Instead, Travel Model geographies could reference census defined areas at the lowest level of aggregation.

Then we just rely on Census to maintain the geometries with no loss of information/accuracy.